% revision from f_run20111211_mixed_sub3items_1fig.m 
% Objective: Customize a visual search model specifically for Scene
% Integratoin


% function [n t]=scene_search(Img, fig_opt)
function f_scene_search(filename, HOMEIMAGES, SalCentroids, SalMask) % note that SalCentroids are x,y index, different from row, column index
% makeSal4SceneInt('D:\Codes\SceneInt\RecogModels\', 'Test_Input.bmp');
% load('sal_output.mat');
v = loadXML(filename); 
VisitedObjs=zeros(1,size(v.annotation.object, 2));

% load image
img_filename = fullfile(HOMEIMAGES, v.annotation.folder, v.annotation.filename); 
img_filename = strrep(img_filename, '/', filesep);
img_filename = strrep(img_filename, '\', filesep);
Orig_Img = imread(img_filename);

flag_subitizing=0;
strategy='mixed';
fig_opt='';
%Img=EyeField(:,:);
confidence_thresh= .9;
noise_sigma=.05;
resp_sigma=.05;
title_text='';
sal_hit_thresh=.6;
sal_terminate_thresh=.4;
%fig_opt=[];
close all;
init;
NLIP_RANGE=30;
Img=reshapeImg(SalMask, field_size); 
SalCentroids=transCentroid(SalCentroids, size(SalMask), field_size);
EyeField=Img;
center=[1+floor(field_size/2), 1+floor(field_size/2)];
center=int16(center);
LIP_Eyepos=zeros(field_size);
LIP_Eyepos(1+floor(field_size/2), 1+floor(field_size/2))=0.6;
total_move=500;
% Layers DLPFC_L2 and DLPFC_L1 compose a recurrent network
% which memorize the location of the items.
DLPFC_L2=zeros(field_size);
DLPFC_L1=zeros(field_size);
DLPFC_L2_fr=zeros(field_size);
DLPFC_L1_fr=zeros(field_size);
DLPFC_inh_mask=zeros(field_size);
DLPFC_inh=zeros(field_size);
DLPFC_Low_THRESHOLD=1.4;
DLPFC_L2_thresh = .1;
LIP_sal=zeros(field_size);
LIP_sal_fr=zeros(field_size);
NLIP=zeros(field_size);
XNLIP=zeros(field_size);
NLIP_inh = zeros(field_size);

% DLPFC Global excitation signal send 
DLPFC_HitFromIT = zeros(field_size);
eta_DLPFC_L2=5;
V2=zeros(field_size);
% confidence_thresh=.5;
% noise_sigma=.01;
% DLPFC_NHi_confidence_sigma;

% removied to fit scene integration
% % NLIP confidence mask
% NLIP_resp_mask=zeros(NLIP_RANGE);
% NLIP_c_idx = int16((NLIP_RANGE+1)/2);
% for i=1:NLIP_RANGE
%     for j=1:NLIP_RANGE
%         if ((i-NLIP_c_idx)^2+(j-NLIP_c_idx)^2)<=((NLIP_RANGE-1)/2)^2
%             NLIP_resp_mask(i,j)=1;
%         end
%     end
% end

% removied to fit scene integration
% % PFC_N(prefrontal) rep of numbers 0-6
size_PFC_N=7;
num_count=0:size_PFC_N-1;
% PFC_NR=zeros(1,size_PFC_N);
% PFC_NS=zeros(1,size_PFC_N);
% PFC_NC=zeros(1,size_PFC_N);
% PFC_NR_fr=zeros(1,size_PFC_N);
% PFC_NS_fr=zeros(1,size_PFC_N);
% PFC_NC_fr=zeros(1,size_PFC_N);
% PFC_NShift=zeros(1,size_PFC_N); % neuron with threshold
% shift_thresh=1;
% shift_amp = 50;
% PFC_NC(1,1)=0.6;
% PFC_NC(1,1)=1;

% for diffenetial equation
dt=.08;
duration_per_move=1.2;
pause_time=.1;

    clim_dlpfc=[0 2];
    clim_vin=[.8 1];
    clim_lip_ctr=[0 2];
    pos=center;mv=[0 0];
    stop_flag=0;  % stop the simulation when stop_flag == 1
    act_record=[];

    %% Initially assume V2<-EyeField
    sim_clk=1;
    if ~isequal(fig_opt, 'fig_off')
        figure; title(title_text); pause(3);
        close all;
        hf_image=figure('position',[620 720  400 300], 'color', 'white');
        h_fig_model=figure('position',[0 0 600 760], 'color', 'white'); title(title_text);
%                         subplot(5,2,2);bar(num_count,PFC_NS);ylabel('DLPFC-N2');
%                         subplot(5,2,4);bar(num_count,PFC_NR);ylabel('DLPFC-N1');
                        subplot(5,2,6);mesh(NLIP);zlabel('LIP-N');
                        subplot(5,2,9);imagesc(1-EyeField, clim_vin);ylabel('Eye Field');colormap('hot');pause(pause_time);
                        subplot(5,2,5); mesh(LIP_sal_fr);ylabel('LIP/FEF Saliency');pause(pause_time);
                        subplot(5,2,7);mesh(LIP_Eyepos);ylabel('LIP Eye Position');  pause(pause_time);  %imagesc(LIP_Eyepos, clim_lip_ctr);
                        subplot(5,2,3); mesh(DLPFC_L1);ylabel('DLPFC-L1');
                        subplot(5,2,1); (DLPFC_L2);ylabel('DLPFC-L2');
%                  h_fig=figure('position',[-1300 -200 600 1200]); %for portait display
    end
    

    %%
    for i=1:total_move
%   dedicate in tracking one object4, avoid distracted by others
%   return inhibition only when target is recognized

        SumLatInh=0;   %DLPFC_N layer 
        %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
         % decide saccade direction vector 'mv'
         % 1. saccade to max of NLIP in subitizing stage
         % 2. saccade to max of LIP_sal in counting stage
         pos=Extract2d_max(LIP_Eyepos);
         

            % Find the current max activation in the eye field
            % * has to be implemented as winner take all later.
%             pos=Extract2d_max(LIP_Eyepos);
            %mv=Extract_NN(LIP_sal_fr, center, sal_terminate_thresh);
            mv=Extract2d_max(LIP_sal_fr);
            % catch the no target exception
            try
                mv=mv-center;
            catch
                mv=[1 1];
            end
            mv=mv./abs(mv);
            mv=int16(mv);
%         end
        %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        
        %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        % at the end of saccade, move eye center away from target to
        % avoid
        % keeping activate LIP neurons
        if stop_flag==1
            %mv=[1 1]; % eye move out the target area
            mv=int16(mv);
        end
        
        % Update eyefield (changed to update by image and current position)
        % to get current position, pos must plus mv
        EyeField=update_eyefield(Img,pos+mv);
        
        if ~isequal(fig_opt, 'fig_off')
            subplot(5,2,9);imagesc(1-EyeField, clim_vin);ylabel('Eye Field');colormap('hot');pause(pause_time);
%             subplot(5,2,2);DrawXML(filename, HOMEIMAGES, );
        end
        %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        
        % check is any target in the visual field,
        % if there is, mark the visited object, and redraw polygons
        if EyeField(center(1), center(2))==1;
            VisitedObjs=Delete_visited(VisitedObjs, SalCentroids, pos+mv);
            RawImg=f_drawXML(filename, HOMEIMAGES, VisitedObjs); 
            figure(hf_image);
            sRawImg=imresize(RawImg.cdata,0.5);
            imshow(sRawImg);
            figure(h_fig_model);
        end 
        
        % Update V2
        V2=EyeField; 
        
        % % compute V2 stimulus to LIP_sal
        % % % this is temporary, since it should either directly saccade to
        % % % destination or get LIP_center info.
        V2_LIP_offset=zeros(field_size);
        V2_LIP_offset(pos(1)+mv(1),pos(2)+mv(2))=1;
        V2_out = conv2(V2, V2_LIP_offset, 'same');
        
        
        % Makes mask that mimics the impulse response of LIP_Eyepos
        w_mask=zeros(3); w_mask(mv(1)+2, mv(2)+2)= 1;
        try
        SubMapLIP = LIP_Eyepos(-1+pos(1):pos(1)+1, -1+pos(2):pos(2)+1);  %the impulse response is very local
        SubMapLIP_fr=mp2fr(SubMapLIP);
        catch
            msg='illegal crop of LIP_Eyepos!'
        end
        
        % Find current peak in DLPFC_N
        % [v, PFC_N_idx] = max(PFC_NHi);

        
        % Dynamic Update
        for j=1:duration_per_move/dt
            
            % %==============================================================
            % compute the local impulse response surrounding the center
            syn_sum=conv2(w_mask, SubMapLIP_fr); syn_sum=syn_sum(2:size(syn_sum,1)-1, 2:size(syn_sum,2)-1);
            Mask_lat_inh_eyepos=ones(3); Mask_lat_inh_eyepos(2,2)=0;
            Lat_inh_sum=conv2(Mask_lat_inh_eyepos, SubMapLIP_fr,'same');
            % Leaky integrator that simulation the dynamic change of the
            % local LIP_Eyepos neurons
            
            SubMapLIP= max(0,SubMapLIP + 5*(-SubMapLIP+syn_sum)*dt);                      %***
            act_record=[act_record max(max(SubMapLIP))];
            % Update the local activation
            LIP_Eyepos(-1+pos(1):pos(1)+1, -1+pos(2):pos(2)+1)=SubMapLIP;
            LIP_Eyepos=max(0,LIP_Eyepos);
            LIP_Eyepos_fr=mp2fr(LIP_Eyepos);
            if ~isequal(fig_opt, 'fig_off')
                subplot(5,2,7);mesh(LIP_Eyepos);ylabel('LIP Eye Position');  pause(pause_time);  %imagesc(LIP_Eyepos, clim_lip_ctr);
            end
            sim_clk=sim_clk+1;  
            %==============================================================
            
            % inhibit subitizing region, when target pos attended (note
            % inhibition taken place only when subitizing is still preceeding, it might disabled when no target found for subitizing ).
            PFC_NS_Input=zeros(1,size_PFC_N);
            max_nlip=max(max(NLIP));
            

                % %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
                % Update activation of DLPFC spatial map
                % Layers DLPFC_L2 and DLPFC_L1 compose a recurrent network
                % which memorize the location of the items.
                % neurons in DLPFC_L2 is not easy to be activated, unless with the
                % cooperating act of LIP_Eyepos and hit(item recognized)
                % signal from IT.

                SynIn = LIP_Eyepos_fr+ceil(LIP_sal(center(1),center(2))-sal_hit_thresh);
                DLPFC_L1 = DLPFC_L1 + dt*(-DLPFC_L1 + 2*DLPFC_L2_fr+max(0,SynIn-DLPFC_Low_THRESHOLD));             %***
                DLPFC_L1_fr=mp2fr(DLPFC_L1);
                DLPFC_L2 = (DLPFC_L2 + dt*(-DLPFC_L2 + DLPFC_L1_fr));                      %***
                DLPFC_L2_fr=mp2fr(DLPFC_L2);

                if ~isequal(fig_opt, 'fig_off')
                   subplot(5,2,3); mesh(DLPFC_L1);ylabel('DLPFC L1'); %imagesc(DLPFC_L1, clim_dlpfc);
                   subplot(5,2,1);mesh(DLPFC_L2);ylabel('DLPFC L2'); %imagesc(DLPFC_L2, clim_dlpfc);
                end

                % Instantiate and update new layer in LIP, which responsible for cobine V2 output
                % and LIP_mem, and by which saccade direction is guided
                DLPFC_inh_mask=zeros(field_size);
                if (DLPFC_HitFromIT(1,1)==1)
                    msg='hit!';
                end
                % revised to increase inhibition range
                % DLPFC_inh_mask(field_size-pos(1)-mv+1,field_size-pos(2)-mv+1)=1;
                DLPFC_inh_mask(field_size-pos(1)-mv(1)-1:field_size-pos(1)-mv(1)+3,field_size-pos(2)-mv(2)-1:field_size-pos(2)-mv(2)+3)=1;
                DLPFC_inh=conv2(max(0,DLPFC_L2_fr-0.05), DLPFC_inh_mask, 'same');
                % %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

            
            % %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            % Updatae Common Activations
            LIP_sal =  max(0, LIP_sal + dt*15*(-LIP_sal + V2 - 3*DLPFC_inh));                    %***
            LIP_sal_fr= mp2fr(LIP_sal);
            
            if ~isequal(fig_opt, 'fig_off')
                subplot(5,2,5); imagesc(LIP_sal_fr);ylabel('LIP/FEF Saliency');
                %subplot(5,2,3); imagesc(DLPFC_inh);pause(pause_time);
            end
        
            % Find the current max activation in the eye field
            % * has to be implemented as winner take all later
            pos_LIP=Extract2d_max(LIP_Eyepos);

            %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

        end 
        
% %     Implement termination saccade mechanism
        tmp=zeros(field_size);
        tmp(pos,pos)=1;
        LIP_sal_tmp=max(0, (LIP_sal_fr - tmp));
        if(max(max(LIP_sal_tmp))<=sal_terminate_thresh)
            if stop_flag==1
                break; 
            else
                stop_flag=1; %stop in next run
            end
        end        
    end
    
    % Report number
%     [v idx]=max(PFC_NR);
%     n=min(6,idx-1);
%     t=i;
end

% Not used, a questionable funtion
function mv=SC(EyeField, center)
    pos=Extract2d_max(EyeField);
    mv =pos-center;
    mv(1)=mv(1)/max(1, abs(mv(1)));
    mv(2)=mv(2)/max(1, abs(mv(2)));    
end

% pp=Eye_field so skiped

function EyeField = update_eyefield(Img,pos)
    % create convolution vector for remapping
    field_size=size(Img,1);
    x=1+field_size-pos(1); % right ward saccade use left ward remapping filter
    y=1+field_size-pos(2);
    mask=zeros(field_size);
    mask(x,y)=1;
    
    % remapping
    EyeField = conv2(Img, mask); 
    l_bound=int16((field_size+1)/2);
    u_bound=field_size+int16((field_size-1)/2);
    EyeField = EyeField(l_bound:u_bound, l_bound:u_bound);
end


function pos=Extract2d_max(M)
    SizeM=size(M);
    Vec=reshape(M, SizeM(1)*SizeM(2), 1);
    [C, idx]=max(Vec);
    pos=[rem(idx(1), SizeM(1)),floor(idx(1)/SizeM(2))+1];
    pos=int16(pos);
end

function target_pos=Extract_NN(Act_map, center, threshold)
    Act_map=Act_map-threshold;
    [row col val]=find(Act_map>0);
    d_row=int16(row)-center(1);
    d_col=int16(col)-center(2);
    Dist=d_row.^2+d_col.^2;
    [v idx]=min(Dist);
    target_pos=[int16(row(idx)) int16(col(idx))];
end

function hold(delay_period)
start_time=clock();
end_time=clock();
    while delay_period>(end_time-start_time)
        end_time=clock();
    end
end

function fr=mp2fr(mp)
fr=1./(1+exp(-7*mp+3.5));
end

function Img_ret=reshapeImg(Img, field_size)
    Img_ret=zeros(field_size);
    Img_ret(field_size/2+1-size(Img,1)/2:field_size/2+size(Img,1)/2,field_size/2+1-size(Img,2)/2:field_size/2+size(Img,2)/2)=Img;
end

function Orig_pos=showFoaImg(pos,center, Orig_Img)
    shift = 18;
    offset=pos-center;
    Orig_center=int16(ceil(size(Orig_Img)/2));
    %FoaImg=Orig_Img(1+(pos(1)-1)*shift:(pos(1)+1)*shift,1+(pos(2)-1)*shift:(pos(2)+1)*shift)
    Orig_pos=Orig_center+offset*shift;
end

function VisitedObjs=Delete_visited(VisitedObjs, SalCentroids, pos)
    for i=1:size(VisitedObjs,2)
        if pos==int16(SalCentroids(i,:))
            VisitedObjs(i)=1;
        end
    end
end

function RetCentroids = transCentroid(SalCentroids, size_sal_msk, field_size)
    shift=int8((field_size-size_sal_msk)/2)
    Shift=ones(size(SalCentroids,1),2);
    Shift(:,1)=shift(1);Shift(:,2)=shift(2);
    RetCentroids=SalCentroids+Shift;
end